/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 Copyright (C) 2008, Proverbio Alessandro and Radice David

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

/**
 * \file solution.cpp
 * \brief This file contains the Solution class
 *
 * The Solution class is to be used to perform post-processing
 * or error estimates on the computed numerical solution
 */

# include "solution.hpp"

namespace spectral {

Solution::Solution(double const * const & u,
                   int element,
                   Mesh const & me) : _weights(NULL),
                                      _base(NULL),
                                      _nnodes(me.getnnodes()),
                                      _size(_nnodes*_nnodes),
                                      _me(me),
                                      _element(element) {
    int k=0;
    int idx;
    _weights=new double[_size];
    _base=new TensorProd[_size];
    for(int i=0;i<_nnodes;++i) {
        for(int j=0;j<_nnodes;++j) {
            idx=_me.getGlobalNum(_element,i,j);
            _base[k]=_me.lagrange(i,j);
            _weights[k++]=u[idx];
        }
    }
}

Solution::Solution(Solution const & u): _nnodes(u._nnodes),
                                        _size(u._size),
                                        _me(u._me) {
    _weights=new double[_size];
    _base=new TensorProd[_size];
    for(int i=0;i<_size;++i) {
        _weights[i]=u._weights[i];
        _base[i]=u._base[i];
    }
}

Solution::~Solution() {
    delete[] _weights;
    delete[] _base;
}

double Solution::operator()(double const & x,double const & y) const {
    double sum=0;
    for(int i=0;i<_size;++i) {
        sum+=_weights[i]*(_base[i])(x,y);
    }
    return sum;
}

void Solution::Dx(int elemnum) {
    for(int i=0;i<_size;++i) {
        _base[i]=_me.Dx(_base[i],elemnum);
    }
}

void Solution::Dy(int elemnum) {
    for(int i=0;i<_size;++i) {
        _base[i]=_me.Dy(_base[i],elemnum);
    }
}

} // namespace
